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Charles W. Therrien and Carlos H. Velasco 
Department of Electrical and Computer Engineering 



Abstract 

A new iterative version of the Proriy method is presented and shown to be 
exceptionally effective in finding ARMA iiio(l('ls for acoustic data in the signal 
domain. The method is based on a quadratic type of gradient algorithm where it 
is shown that the gradient and Hessian are easily computed from the data. The 
new algorithm is found experimentally to have excellent convergence behavior. 
The performance of the algorithm is demonstrated and compared to that of the 
basic Prony method and to that of the Steiglitz and McBride iterative prefiltering 
algorithm on some recorded acoustic data. 



1 Introduction 

Linear pole-zero (ARMA) models are used for a variety of applications in signal pro- 
cessing. These applications include speech and image coding, spectrum estimation, and 
others. Both deterministic and stochastic methods for signal modeling have been stud- 
ied e.xtensively in the literature. Typically deterministic models attempt to represent 
the signal as the impulse response of a suitably chosen linear system while stochastic 
methods attempt to model the signal as the response of the linear system to white noise. 
However in many cases the procedures embedded in these methods (such as estimating 
a correlation matrix for the data) are similar so the philosophical distinction between 



“This work weis sponsored by the Office of Naval Research. 



I 



deterministic and stochastic models is of less significance. A review of various methods 
can be found in many places such bls [l, 2, 3]. 

In our apphcation we are interested in reproducing the time domain w’aveform as 
accurately as possible (as meeisured e.g., by the sum of squared errors). To do this 
frequently requires modeling procedures that can produce nonminimum-phase transfer 
functions [4, 5] and this eliminates most of the stochastic methods. Consequently, it is 
most appropriate here to think of the signal as being modeled as the impulse response 
of a linear system and seek to minimize the sum of squared errors between this impulse 
response and the given data. 

A well known approach to the modeling problem is Prony’s method, w'here the data 
is represented in the time domain b}’ a sum of damped exponentials w'ith appropriate 
weighting coefBents. Prony’s method avoids the direct nonlinear problem of minimizing 
the sum of squared errors and solves a pair of linear problems. However in many cases 
involving real data this results in a suboptimal model as shown for the acoustic data 
(corresponding to the ringing of a wrench hit on the floor) in Fig. 1(a). In this paper we 
present a method for iteratively adjusting the parameters (poles and zeros) in Prony’s 
method to values that further minimize the sum of squared errors. The method, which 
we call the iterative Prony method, produces a much improved model cls shown in Fig. 

1(b). 

Iterative methods for ARMA modeling are not new’ [6, 7. 8, 9. 10]. Most stochastic 
methods inv’olving maximum likelihood solutions involve iteration (see e.g., [6]) and a 
very effectiv’e method known as iterative prefiltering [7, 8] is capable of results similar 
to the iterative Prony method. The latter, like many other methods (e.g., [9]), uses the 
polynomial coefficients ais the parameters and thus does not have a clear interpretation 
in terms of placement of poles and zeros. Iterative prefiltering can also lead to oscillatory 
behavior in the sum of squared errors while the iterative Prony method tj'pically results 
in a monotonic decrease of the error norm and seems to have better behavior when the 
model order is not knowm. 

The remainder of this paper is organized as follows. Section 2 states the basic form of 
the modern Pronj’ method and estabhshes our notation. Section 3 develops the iterative 
Prony algorithm. Section 4 gives an example and comparison of performance on recorded 
acoustic data. Section 5 provides a summary and conclusions. 
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Figure 1: Approximateh' 10 ms segment of the sound corresponding to the ringing of 
a wrench hit on the floor, (a) Model produced by the (non-iterative) Prony method (4 
poles and 3 zeros), (b) Model produced by the iterative Pron}' method (also 4 poles and 
3 zeros). 
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2 Prony’s Method of ARMA Modeling 

Although there are many variations (see [11, Chapter 11]). the modern Prony method 
can be thought of as a method for approximating a data sequence x[n] by a sequence 
x[n] which is the impulse response of a rational hnear system satisfying the difference 
equation 

x[n] + ajxfn — 1] + • • • + apx[n — P] = + biS[n — l] + • • • + bgSln — Q\ 

( 1 ) 

If the requirement x[n] = x[nj; n = 0, 1, A’j — 1 is applied to (1), where is the 

length of the data, and the difference equation is evaluated for n = 0. 1, . . . , — 1, the 

result is the matrix equation 
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which can be written in partitioned form as 




Thus if A’i > P + (5 + 1 a least squares solution for the a, ( AR) parameters can be found 
from the lower set of partitions, and the upper set of partitions can be used to find the 
bj (MA) parameters. This is the fundamental idea underlying Prony 's method. 

Prony’s method is frequently expressed directly in the signal domain. Instead of 
solving (3) for the vector b, we can instead find the roots of the denomimator polynomial 
of the system in (1) 

A{z) = \ + + 02 Z~‘^ -\ \- apz~^ (4) 
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and express the impulse response i[n] in terms of these roots. Specifically if Tj, rj, . . . , rp 
are the roots of (4) (assumed to be distinct) then x[n] can be written as 



x[n] = cir" + C 2 rj 4 h CpVp 



(5) 



Again, if we require i[n] = i[n]; n = 0,1 Ns — I, then by evaluating (5) for 

n = 0, 1, . . . , A’^s — 1 we have the matrix equation 
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( 6 ) 



which can be solved in a least squares sense to obtain the coefficients c,. This implemen- 
tation of Prony's method has an advantage over simply solving for b from (3) since all 
of the data is used in the computation. When (2) or (3) is used to find the bj, only the 
data values from 0 to Q axe used in the computation. 

In the case of multiple roots at the same location a slight variation of the same 
procedure can be used. Suppose, for example, that rj is a double root. In this case, x[n] 
has the form 

x[n] = Cir" + c^nrj -f • • • -1- cpVp (7) 

and the matrix equation to solve for the coefficients becomes 
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( 8 ) 



This situation is rare however, because computational errors and errors inherent to the 
modeling method itself contribute to produce roots that may be very close to each other 
but not exactly at the same location. 

While Prony’s method represents a clever scheme to separate the solution of the AR 
and MA parameters in the modeling problem, the separation is frequently achieved at the 
cost of a degradation in performance, as mentioned in the introduction. It is well known, 
for example, that Prony’s method does not choose the AR parameters to minimize any 
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norm of the error e[n] = x[n] — x[n] but rather chooses them to minimize the norm 
of a modified error of linear prediction resulting from processing the original data x[n] 
through the filter represented by A{z). As a result, there is motivation to modify the 
model obtained by Prony’s method through an iterative scheme that seeks to directly 
minimize the norm of the error e[n]. This method retains the separabifity of the original 
Prony method and so is called the iterative Prony method. It differs from maximum 
likelihood in that there is no stochastic criterion and it does not attempt to solve for the 
AR and MA parameters simultaneously. The method is developed below. 

3 An Iterative Prony Algorithm 

3.1 Iterative Algorithms 

A multidimensional function Q(ri, rj, • • • , rp) that is continuous and differentiable can be 
minimized using one of several very powerful optimization techniques known as gradient 
methods [12]. Some of those methods are derived on the basis of a quadratic model that 
can be obtained from a truncated Taylor series expansion of Q(r). Let denote the 
value of r at the iteration. Then for any point r = + 6 when 6 is small, the 

function can be approximated by 

where g and G represent the vector of first derivatives (gradient) and the matrix of second 
derivatives (Hessian) of the function Q(r) and they should be available at every point. 
In Newton’s method the iterate is taken to be where the correction 

minimizes q^^^{6). This method is only well defined when the Hessian matrix G 
is positive definite, in which case the iteration of Newton’s method is given by the 
following procedure [13]: 

1. solve = -gW for 6=6^'^^ 

2. set + 6 

The fact that G^^l may not be positive definite when x^*^ is far from the solution, and 
that even when G*^'^ is positive definite convergence may not occur, makes this method 
undesirable as a general formulation of a minimization algorithm. 

Since Newton's method is defined only when the matrix G**^* is positive definite, and 
this matrix is positive definite only when the error 6 is “small”, we can say that the 



6 



method is defined only in some neighborhood of in which ) agrees with 

Q(r(*=) -f 6) in some sense. In such cases, it is correct to choose with 

the correction <5*^^ minimizing for all ^ in fi**^*. This method is referred to 

35 the restricted step method because the step is restricted by the region of validity of the 
Taylor series [13]. The region of definition for the iteration can then be expressed a5 

= {r : llr-r'^^^ll < h<^>} (11) 

where || • || denotes the norm of the vector. In this case, the optimization problem can 
be stated as: 

minimize^ q^^\s) subject to ||6|| < (12) 

As mentioned before, the least squares norm || • jjj is the one most commonly used in 
these type of problems so it is the one used in this paper. The problem that now becomes 
apparent is how to select the error margin of the neighborhood (11). This margin 
should be as large as possible because the iteration step is directly related to it. Various 
methods have been proposed to control the parameter One of these methods [13] 

attempts to insure that the Newton's minimization problem (10) is always defined b}’ 
adding a multiple of the unit matrix I to and computing the new problem 

+ (13) 

where the net effect is that increases in v cause ||6||2 to decrease, and vice v'ersa. 

If we define 

of*') = — ! L (14) 

^ A^f'') Q{k] - V ) 

then the ratio represents a measure of the accuracy to which approximates 

0(rf'') + (5^*^^) on the k''^ step, and as the accuracy increases pf'^) gets closer to unity. 
Using (14), Marquardt [14] suggested an algorithm that tries to adaptively maintain 
as large as possible while controUing the ratio p**'). The k^^ iteration of such an algorithm 
is stated as 

1. given and calculate gf^) and G^*^); 

2. factor G^*^) + if not positive definite, reset and repeat; 

3. solve (13) to find 

4. evaluate + 6**^)) and hence pf'’'); 



5. if < 0.25 set = 4i/W 
else if > 0.75 set 

else set 

6. if < 0 set else set 

Here the parameters 0.25, 0.75, 4 and 2 are arbitrary and > 0 is also chosen arbi- 
trarily. Proofs of global and second order convergence for this algorithm are given in [13] 
for the cases when the first and second derivatives of the function Q(r) exist, and the 
vector belongs to a bounded n-dimensional space for all k. Although this method 
does have some disadvantages, it represents a good basis for the formulation of a general 
minimization algorithm. 



3.2 Computation of Derivatives 

Let us now return to the problem of representing a sequence x[n] as a linear combination 
of complex exponentials. Equation 6 can be written as 



Rc = X -f c 



(15) 



where e is the equation error, x represents the data, which may or may not be complex, 
c is the vector of complex coefficients, and R is the matrix of complex roots, which can 
be written in terms of real and imaginary parts as 
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By defining 



Q(r) = ||e ||2 = = (Rc - x)*^ (Rc - x) 



(17) 



it is clear that the problem is to find the vector r = r/j 4- jr; of P complex roots that 
minimizes the function Q{r). 
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To make the following development less cumbersome let us define the gradient oper- 
ator with respect to the real and imaginary parts of the roots as 

d 
8 



Vr = 
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V, 



rR 
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drj. 
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(18) 



Then the gradient vector g of first derivatives and the Hessian matrix G of second 
derivatives can be written compactly as 



g = VrQ = 



Vr,C 



and 



G = V, [(V,2)’ 






T 1 



(19) 



( 20 ) 



Vr,(V,„0)' V„(V„S)' 

Note that it is essential to work with the derivatives involving the real and imaginary 
parts of the roots rather than the complex roots themselves because Q(T) in (17) is not 
an analytic function of the complex variables rj , T 2 , . . . . rp. It is shown in the appendix 
that the upper and lower partitions of g are given by 

Vr„C = 2Re{F-’-«) (21) 

and 

= 2 Im{F'^} (22) 

where F is a matrix with columns f,- defined by 
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Further, if S is defined as the matrix with columns 
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(24) 



then it can be shown (see appendix) that the partitions of G are given by 



= 2Re {F*"F + diag (S*"e)} (25) 

= 2Im {F'^F - diag (S*^)} (26) 

Vr;(Vr«Q)^ = -2Im {F*"F - diag (S*"^)} (27) 

Vr,(Vr,Q)^ = 2Re {F'^F - diag {S'^e)} (28) 



where the notation “diag ” applied to a vector represents the operation of forming a 
diagonal matrix w'hose components are the components of the vector. Thus the gradient 
g and Hessian G are convenient to compute. 



3.3 Real Axis Poles and Zeros 

Some discussion is necessary about how the iterative Prony algorithm handles poles and 
zeros on the reaJ axis. \\'hen such poles and zeros occur, the imaginary parts and all of 
the associated derivatives are zero. Thus poles and zeros initiall}' on the real axis remain 
on the real axis through all successive iterations. This can be a practical disadvantage 
because the original Prony algorithm, which is used to produce the starting configuration 
of poles and zeros, can sometimes place a pair of poles or zeros on the real axis, while in 
fact this pair really belongs off of the real axis in a conjugate symmetric configuration. 
In this situation the method so far described will never reach the true configuration. A 
modification is therefore introduced to defiberately displace those roots from the real 
axis. The algorithm places the roots in a conjugate symmetric position by calculating 
the mean of their values and adding to them a small arbitrary offset in the imaginary 
direction, and proceeds with the iterations. K the tendency of the roots is to “go back” to 
the real axis then they are returned to their initial position and the iterations continue. 
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Figure 2: Displacement of the poles and zeros of an order (4,3) model; the modeled signal 
in this case actually heis two poles on the real axis. 

Otherwise the roots ma\' continue to spread apart and move as a complex pair. Figure 
2 is an example of the displacement of the poles and zeros of an order [P.Q) — (4,3) 
model. In this case the modeled signal actually has two poles on the real axis and the 
initial model correctly placed two of the poles on the real axis. Those poles are displaced 
from the real axis by an amount of approximately 0.3 in the imaginary direction by the 
algorithm, but then after some iterations it is clear that the poles are tending to return 
to the real axis. At this point the poles are returned to the real axis by setting their 
imaginary parts to zero and the iterations continue until convergence is obtained. The 
opposite situation is shown in Figure 3. In this case the initial model also has two poles 
located on the real axis, but contrary to the case presented above, the roots, after being 
displaced from the real axis, (again by an amount of 0.3), continue to move away from 
the cLxis until they reach their final position in the first and fourth quadrants closer to 
the unit circle. 
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Figure 3: Displacement of the poles and zeros of an order (4.3) model. The initial model 
shown has poles on the real axis: the final model in this case has those poles moved to 
conjugate symmetric locations. 
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4 Performance Comparison 

The iterative Prony method has been tested on a large variety of data including both 
simulated and recorded sonar and speech data. In testing, we compared the performance 
of the new iterative Prony algorithm to that of the iterative prefiltering algorithm, which 
previously had been found to be the algorithm most suitable for our work in waveform 
modehng of acoustic data in the signal domain [5]. In tests on simulated data, we found 
that when the correct order of the model is known and it is used to model the data, then 
both iterative algorithms tested produced very good models. However, when given an 
order higher than the true order of the signal, the iterative prefiltering algorithm tended 
to oscillate and took longer to reach convergence. On the other hand, in most cases the 
performance of the iterative Prony method was not affected by this change in the order 
of the model. The final model w'ould simply have extra poles and zeros positioned to 
cancel each other. 

On the real recorded data we found that in all cases the iterative Prony method 
provides much better models than the regular (non-iterative) Prony method. In many 
cases, the performance of the iterative prefiltering method was comparable to that of 
the iterative Prony method, but in a large number of cases the error produced by the 
iterative prefiltering algorithm was highly oscillatory. The same number of iterations 
was used for both the iterative Prony and the iterative prefiltering methods in all of the 
tests. In order to obtain the best possible results from the iterative prefiltering algorithm 
(especially in those cases when there was no convergence), the model selected was that 
corresponding to the iteration with the lowest error between the model and the original 
sequence. 

Figure 4 shows the results of modeling a 100-point segment (approximately 8 ms) 
of underwater recorded data (ice cracking). The signal was modeled using the signal 
domain Prony method of section 2 and the two iterative algorithms mentioned above. A 
model of order (12,11) was used in all cases. It is clear that the models produced by the 
iterative methods (Fig. 4(b) and (c)) represent the original signal much more faithfully 
than the model produced by the non-iterative method (Fig 4(a)). It can also be seen that 
the model produced by the iterative Prony method provides a much better representation 
of the original signal than the model produced using the iterative prefiltering algorithm. 
Figure 5 shows the behavior of the sum of squared errors at each iteration between the 
original signal and the models produced by the two iterative methods compared here. 
The oscillating behavior of the error when using the iterative prefiltering algorithm was 
found to be one of the main disadvantages of using this method. Because the iterative 
Prony method is based in minimizing the error between the model and the original signal 
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magnitude magnitude magnitude 





Figure 4: Segment of recorded underwater sound corresponding to ice cracking and three 
ARMA models, (a) Non-iterative Prony model, (b) Iterative prefiltering model, (c) 
Iterative Prony model. 
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Figure 5: Sum of squared errors at each iteration between the original and the models 
produced by the iterative prefiltering and iterative Prony methods. 

by directly adjusting critical parameters, the behavnor of the error tends to be smoother 
and to decrease monotonically toward convergence with increasing number of iterations. 
In this case the error rises slightly at the third iteration but then decreases monotonically 
for the remainder of the test. 



5 Conclusions 

ARMA models can produce very accurate reproductions of short segments of acoustic 
data but the accuracy depends critically on the method used to find the model parame- 
ters. Frequent!}' nonminimum-phase models are required to match the data in the signal 
domain. Prony’s method is convenient to use for ARMA modehng because it finds the AR 
and MA parameters separately by solving linear equations. However the model provided 
by Prony’s method is frequently suboptimal. The iterative Prony method described in 
this paper results in a significant performance improvement over the basic Prony method 
in practical problems while retaining the separability characterizing the original method. 
The new method is based on iteratively moving the poles of the model in directions that 
tend to decrease the error between the original and modeled data and has been found 
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experimentally to have excellent convergence properties. 

The new method ha^ also been compared with the iterative prefiltering algorithm of 
Steiglitz and McBride [7]. Neither method is guaranteed to have monotonic convergence 
of the error, but in many cases of testing on real and simulated data, where the true model 
order was not known, the iterative prefiltering algorithm had an oscillatory behavior 
while the iterative Prony method showed monotonic or near monotonic convergence. The 
computational requirements of the two algorithms were compared by testing and counting 
floating point operations within the MATLAB implementation. For a complex data set, 
the computational requirements for the iterative Prony method are approximately 

672P^ + (24A^, + 102)F^ + (60A^, + 46)P + 198A', 

where P is the number of poles and Ns is the length of the data, while those for the 
iterative prefiltering algorithm are 

64(P + Q - If + SNs{P + Q? + 10(P + Q)Ns + 12AT 

where Q is the number of zeros. In typical problems of models in the range of 8 to 16 
poles and zeros the iterative Prony method requires about 25% more operations than 
the iterative prefiltering method. However the likelihood of the iterative Prony method 
to reach convergence in practical applications can result in running the algorithm for a 
fewer total number of iterations and thus an overall decrease in the total computation. 
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A Appendix 



A.l Gradient vector g 

From (15) it is easy to see that 



d 



dR 



c and 

dr, dr, 



de’-^ 

dr. 



= c 



• T 



dr. 



(29) 



where r, here represents the real or imaginary part of the root. Using (29) and the 
chain rule in (17) leads to 



dQ 

drR^ 



= € 



■ gR 

drR^ 



c + c 



.r dR’^ 

drR, 



€ = 2Re 



- aR 

drn^ 



(30) 



and explicit evaluation of the partial derivative of the matrix R results in 



dR 



dr 



-c = 



R. 



0 
1 

2(r«, + j>;.) 

3(r^, +jri,f 

. {Ns - 1) (rR, 



def f 

C, = f, 



(31) 



where c, is the component of the vector c. Then, using (31). equation 30 can be 
written as 

do 

^ = 2Re = 2Re (32) 

drR^ 

Finally, using (32) for i = 1 • • • P, the upper partition of (19) becomes 



(33) 







r fr 1 




Vr^Q = 2Re< 




fr 


€ ^ 






1 


J 



which is the same as (21). 

A similar procedure can be used to obtain the gradient of Q with respect to the 
imaginary part of the vector r. Once more, from (29) and the chain rule applied to (17) 
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it follows that 



Thus 



dQ 

drj, 



n 

II 


(34) 


ori. 


(35) 


j€-^f,-jfr€ = 2Im [f^e] 


(36) 



and (22) follows. 

A. 2 Hessian matrix G 

An expression for the Hessian matrix G can be obtained as follows. Using a somewhat 
more convenient notation (20) can be rewritten as 



d 


r ^ 

[ SrK, 


^ ■ 




— 1 


■ ^ 
C"’H, 


!£. 1 





z = L...,P 
k = l,...,P 



G = 



then substituting (32) and (36) into (37) yields 



G = 



Now using the chain rule and both expressions in (29) leads to 



(37) 



<^TR^ \ 


' + f-'e] 


; - f.-'e| ) ■ 


z = l,.. 


..P 


1 

a^| 

pr 




i[«-f,-f-.|) _ 


Ar = l,.. 


.,P. 



(38) 



G = 



+ f; 



• T ^ I 



T C7XV I 

^ drji^ * ‘t dTji^ ~ ’ 



T df, I _ r*T dK p _ 

dTR^ dvR^ S drR^ dn 






*T dii 
dri. 



9rR. drR 



+ C 



»T r ^R df*’^ 



z = 1, . . . , P; A: = 1, . . . , P 



(39) 
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and from the definition of f, in (31) it is seen that 

0 
0 



df, 






6(r/i, +;>/.) 

12 (rR. +jri,f 

(A; -1)(A^- 2) (rR. +;>;.) 



N.-3 



c, if i=k 



(40) 



0 



otherwise 



= s,6,k 



Thus if s, is defined by (24) it is seen that 
where 6,* is the dirac delta. In the same way it can be shown that 

dU 



dru 



= js^Sa 



(41) 



(42) 



These last two expressions together with (31), (34). and (35) can now be substituted in 
(39) to obtain 

■ j - frf* - s^eS,,] • 

G = 

_ j [e-^s.<5., - - f-f, + s‘JeS,,] . 

f = l,...,P; ^:=1,....P. (43) 

Since the elements of the matrix G are formed by additions and subtractions of complex 
scalars with their respective complex conjugates, an alternate expression for G is 

• 2Re + 2Re 2Im - 2Im [s'JeS,,] 

G = 

_ -2Im + 2Im [5^6 4*] 2Re - 2Re [sre6.-,] 

f = 1,...,P; k = \ (44) 

Finally, since the notation in (44) specifies the element of each partition, we notice 
that it is equivalent to (25) through (28). 
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